In this notebook, we try to visualize the automatic occupancy detection algorithm used in the RMV2.0 time-of-week temperature (TOWT) model. This article does not promote use of the algorithm for all cases. You should determine whether it is appropriate for your building and usage history.
The occupancy detection is done by findOccUnocc(). Comments in the code explain the proccess:
Define ‘occupied’ and ‘unoccupied’ based on a regression of load on outdoor temperature: times of week that the regression usually underpredicts the load will be called ‘occupied’, the rest are ‘unoccupied’ This is not foolproof but usually works well.
If the regression underpredicts the load more than 65% of the time then assume it’s an occupied period.
Some details are important to clarify:
Here is the call stack:
train_model (called at server.R#298, defined at utils.R#215)towt_baseline (called at utils.R#243, defined at towt_baseline.R#42)makeBaseline (called at towt_baseline.R#74, defined at towt_baseline.R#433)fitLBNLregress (called at towt_baseline.R#491, defined at towt_baseline.R#224)findOccUnocc (called at towt_baseline.R#296, defined at towt_baseline.R#186)Can we recover enough data to show the process, from a saved project file? If not, at what point in the code do we need to store some more data?
library(dplyr)
library(glue)
library(RMV2.0)
# The 'project' file saved when you used the RMV2.0 add-in.
rds_file <- "C:/RMV2.0 Workshop/something/Project_02.19.rds"
Project <- readRDS(rds_file)
i = 1
What next? Here, we recreate the function call to towt_baseline.
pre_Data_i <- Project$Data_pre[[i]]
timescaleDays <- Project$model_obj_list$models_list[[i]]$towt_model$timescaleDays
intervalMinutes <- Project$Data_pre_summary[1,6]
fahrenheit <- Project$fahrenheit
res_baseline <- towt_baseline(train_Data = pre_Data_i,
pred_Data = pre_Data_i,
timescaleDays = timescaleDays,
intervalMinutes = intervalMinutes,
fahrenheit = fahrenheit,
)
Loading required package: shiny
For time weights, tcenter = 2006-01-01 01:00:00 MST
In fitLBNLregress, reduced tempKnots:
[1] 4.444444 12.777778 18.333333 26.666667
For time weights, tcenter = 2006-03-15 01:00:00 MST
In fitLBNLregress, reduced tempKnots:
[1] 4.444444 12.777778 18.333333 26.666667
For time weights, tcenter = 2006-05-27 MST
In fitLBNLregress, reduced tempKnots:
[1] 4.444444 12.777778 18.333333 26.666667
For time weights, tcenter = 2006-08-08 MST
In fitLBNLregress, reduced tempKnots:
[1] 4.444444 12.777778 18.333333 26.666667
For time weights, tcenter = 2006-10-19 23:00:00 MST
In fitLBNLregress, reduced tempKnots:
[1] 4.444444 12.777778 18.333333 26.666667
For time weights, tcenter = 2006-12-31 23:00:00 MST
In fitLBNLregress, reduced tempKnots:
[1] 4.444444 12.777778 18.333333 26.666667
Now, to illustrate some of the data structures used in the next step …
train <- Project$model_obj_list$models_list[[i]]$train
train
train$time <- as.POSIXlt(train$time, format="%m/%d/%y %H:%M")
head(train$time)
[1] "2006-01-01 01:00:00 MST" "2006-01-01 02:00:00 MST" "2006-01-01 03:00:00 MST"
[4] "2006-01-01 04:00:00 MST" "2006-01-01 05:00:00 MST" "2006-01-01 06:00:00 MST"
pred <- Project$model_obj_list$models_list[[i]]$pred
pred
pred$time <- as.POSIXlt(pred$time, format="%m/%d/%y %H:%M")
head(pred$time)
[1] "2006-01-01 01:00:00 MST" "2006-01-01 02:00:00 MST" "2006-01-01 03:00:00 MST"
[4] "2006-01-01 04:00:00 MST" "2006-01-01 05:00:00 MST" "2006-01-01 06:00:00 MST"
What next? Here, we recreate the function call to makeBaseline. The function loops over a list of timestamps spaced by timescaleDays (the hyperparameter set by the user in the GUI, from 15 to 90 days). At each step, it sets up weights centered at the timestamp and calls fitLBNLregress to create a mini-regression. It then bundles all these regressions together.
verbosity=5
towt_model <- makeBaseline(train$time,
train$eload,
train$Temp,
pred$time,
pred$Temp,
intervalMinutes=intervalMinutes,
timescaleDays=timescaleDays,
fahrenheit=fahrenheit,
verbose=verbosity)
[1] "starting makeBaseline()"
[1] "running regression at 6 steps"
[1] "starting model run number 1"
For time weights, tcenter = 2006-01-01 01:00:00 MST
[1] "starting fitLBNLregress()"
[1] "done determining occupied hours"
In fitLBNLregress, reduced tempKnots:
[1] 4.444444 12.777778 18.333333 26.666667
[1] "done setting up temperature matrix"
[1] "fitting regression for occupied periods"
[1] "done with prediction for occupied periods"
[1] "fitting regression for unoccupied periods"
[1] "done with prediction for unoccupied periods"
[1] "leaving fitLBNLregress()"
[1] "starting model run number 2"
For time weights, tcenter = 2006-03-15 01:00:00 MST
[1] "starting fitLBNLregress()"
[1] "done determining occupied hours"
In fitLBNLregress, reduced tempKnots:
[1] 4.444444 12.777778 18.333333 26.666667
[1] "done setting up temperature matrix"
[1] "fitting regression for occupied periods"
[1] "done with prediction for occupied periods"
[1] "fitting regression for unoccupied periods"
[1] "done with prediction for unoccupied periods"
[1] "leaving fitLBNLregress()"
[1] "starting model run number 3"
For time weights, tcenter = 2006-05-27 MST
[1] "starting fitLBNLregress()"
[1] "done determining occupied hours"
In fitLBNLregress, reduced tempKnots:
[1] 4.444444 12.777778 18.333333 26.666667
[1] "done setting up temperature matrix"
[1] "fitting regression for occupied periods"
[1] "done with prediction for occupied periods"
[1] "fitting regression for unoccupied periods"
[1] "done with prediction for unoccupied periods"
[1] "leaving fitLBNLregress()"
[1] "starting model run number 4"
For time weights, tcenter = 2006-08-08 MST
[1] "starting fitLBNLregress()"
[1] "done determining occupied hours"
In fitLBNLregress, reduced tempKnots:
[1] 4.444444 12.777778 18.333333 26.666667
[1] "done setting up temperature matrix"
[1] "fitting regression for occupied periods"
[1] "done with prediction for occupied periods"
[1] "fitting regression for unoccupied periods"
[1] "done with prediction for unoccupied periods"
[1] "leaving fitLBNLregress()"
[1] "starting model run number 5"
For time weights, tcenter = 2006-10-19 23:00:00 MST
[1] "starting fitLBNLregress()"
[1] "done determining occupied hours"
In fitLBNLregress, reduced tempKnots:
[1] 4.444444 12.777778 18.333333 26.666667
[1] "done setting up temperature matrix"
[1] "fitting regression for occupied periods"
[1] "done with prediction for occupied periods"
[1] "fitting regression for unoccupied periods"
[1] "done with prediction for unoccupied periods"
[1] "leaving fitLBNLregress()"
[1] "starting model run number 6"
For time weights, tcenter = 2006-12-31 23:00:00 MST
[1] "starting fitLBNLregress()"
[1] "done determining occupied hours"
In fitLBNLregress, reduced tempKnots:
[1] 4.444444 12.777778 18.333333 26.666667
[1] "done setting up temperature matrix"
[1] "fitting regression for occupied periods"
[1] "done with prediction for occupied periods"
[1] "fitting regression for unoccupied periods"
[1] "done with prediction for unoccupied periods"
[1] "leaving fitLBNLregress()"
[1] "leaving makeBaseline()"
What next? Here, we recreate the function call to fitLBNLregress.
dataTime <- train$time
dataLoad <- train$eload
dataTemp <- train$Temp
predTime <- pred$time
predTemp <- pred$Temp
tempKnots = (c(40, 55, 65, 80, 90)-32)*5/9
doTemperatureModel<-T
verbose<-5
npoints = length(dataLoad)
t0 = min(dataTime,na.rm=T)
t1 = max(dataTime,na.rm=T)
deltaT = as.numeric(difftime(t1,t0,units="days"))
nsegments = max(1,ceiling(deltaT/timescaleDays))
segmentwidth = (npoints-1)/nsegments
pointlist = floor(sort(npoints-segmentwidth*(0:nsegments))+0.001)
nModelRuns = max(1,length(pointlist))
#for (irun in 1:nModelRuns)
irun <- 1
tcenter = dataTime[pointlist[irun]]
tDiff = as.numeric(difftime(tcenter,dataTime,units="days"))
tDiffPred = as.numeric(difftime(tcenter,predTime,units="days"))
weightvec = timescaleDays^2/(timescaleDays^2 + tDiff^2)
regOut = fitLBNLregress(dataTime, dataLoad, dataTemp, predTime, predTemp,
tempKnots = tempKnots, weightvec=weightvec,
intervalMinutes=intervalMinutes,fahrenheit=fahrenheit,
doTemperatureModel=doTemperatureModel,verbose=verbose)
[1] "starting fitLBNLregress()"
[1] "done determining occupied hours"
In fitLBNLregress, reduced tempKnots:
[1] 4.444444 12.777778 18.333333 26.666667
[1] "done setting up temperature matrix"
[1] "fitting regression for occupied periods"
[1] "done with prediction for occupied periods"
[1] "fitting regression for unoccupied periods"
[1] "done with prediction for unoccupied periods"
[1] "leaving fitLBNLregress()"
What next? Here, we recreate the function call to findOccUnocc (from the first step inside fitLBNLregress). Note:
intervalOfWeek using the weekday, hour, and minute parts of the timestamp data.intervalOfWeek increments by 1 every intervalMinutes (in this project, every 60 minutes).intervalOfWeek = 1.timeVec <- dataTime
loadVec <- dataLoad
tempVec <- dataTemp
#predTime
#predTemp
#tempKnots
#weightvec
#intervalMinutes
#fahrenheit
#doTemperatureModel
#verbose
minuteOfWeek = 24*60*timeVec$wday+60*timeVec$hour + timeVec$min
intervalOfWeek = 1+floor(minuteOfWeek/intervalMinutes)
# If we have temperature data then fit the time-of-week-and-temperature model
if (fahrenheit) {
# temperature vector is already in fahrenheit
tempVecF = tempVec
tempVec = (tempVec-32)*5/9
tempVecPredF = predTemp
tempVecPred = (predTemp-32)*5/9
} else {
tempVecF = (tempVec*9/5)+32
tempVecPredF = (predTemp*9/5)+32
tempVecPred = predTemp
}
# findOccUnocc requires Fahrenheit temperatures; everywhere else we can use either
# Celsius or Fahrenheit, as long as temperature knots are set appropriately
#
# base occupied/unoccupied decision only on cases where we have load data:
okload = !is.na(loadVec)
occInfo = findOccUnocc(intervalOfWeek[okload],loadVec[okload],tempVecF[okload])
head(occInfo,40)
uTOW okocc
[1,] 2 0
[2,] 3 0
[3,] 4 0
[4,] 5 0
[5,] 6 0
[6,] 7 0
[7,] 8 0
[8,] 9 0
[9,] 10 0
[10,] 11 0
[11,] 12 0
[12,] 13 0
[13,] 14 0
[14,] 15 0
[15,] 16 0
[16,] 17 0
[17,] 18 0
[18,] 19 0
[19,] 20 0
[20,] 21 0
[21,] 22 0
[22,] 23 0
[23,] 24 0
[24,] 25 0
[25,] 26 0
[26,] 27 0
[27,] 28 0
[28,] 29 0
[29,] 30 0
[30,] 31 0
[31,] 32 0
[32,] 33 1
[33,] 34 1
[34,] 35 1
[35,] 36 1
[36,] 37 1
[37,] 38 1
[38,] 39 1
[39,] 40 1
[40,] 41 1
tail(occInfo,20)
uTOW okocc
[149,] 150 0
[150,] 151 0
[151,] 152 0
[152,] 153 1
[153,] 154 1
[154,] 155 1
[155,] 156 1
[156,] 157 1
[157,] 158 1
[158,] 159 0
[159,] 160 0
[160,] 161 0
[161,] 162 0
[162,] 163 0
[163,] 164 0
[164,] 165 0
[165,] 166 0
[166,] 167 0
[167,] 168 0
[168,] 1 0
What next? Here, we demonstrate the occupancy detection algorithm within findOccUnocc. Note:
uTOW = unique time-of-week. If intervalMinutes=60, then there are up to 168 such numbers. We do not need to sort them, which is unnecessary.intervalOfWeek2 <- intervalOfWeek[okload]
loadVec2 <- loadVec[okload]
TempF <- tempVecF[okload]
#intervalMinutes
#verbose
# Figure out which times of week a building is in one of two modes
# (called 'occupied' or 'unoccupied')
uTOW = unique(intervalOfWeek2)
nTOW = length(uTOW)
# Define 'occupied' and 'unoccupied' based on a regression
# of load on outdoor temperature: times of week that the regression usually
# underpredicts the load will be called 'occupied', the rest are 'unoccupied'
# This is not foolproof but usually works well.
#
TempF50 = TempF-50
TempF50[TempF > 50] = 0
TempF65 = TempF-65
TempF65[TempF < 65] = 0
if (verbose > 4) {
cat("Fitting temperature regression...\n")
}
Fitting temperature regression...
amod = lm(loadVec2 ~ TempF50+TempF65,na.action=na.exclude)
okocc = rep(0,nTOW)
cat("Detecting occupancy ...\n")
Detecting occupancy ...
for (itow in 1:nTOW) {
okTOW = intervalOfWeek2==uTOW[itow]
# if the regression underpredicts the load more than 65% of the time
# then assume it's an occupied period
if ( sum(residuals(amod)[okTOW]>0,na.rm=T) > 0.65*sum(okTOW) ) {
okocc[itow]=1
}
cat(glue('[{format(t,width=3)}] {format(nunder,width=2)}/{format(ntotal,width=2)} -> {occ}',
t=uTOW[itow],
nunder=sum(residuals(amod)[okTOW]>0,na.rm=T),
ntotal=sum(okTOW),
occ=okocc[itow]
),'\n')
}
[ 2] 0/53 -> 0
[ 3] 0/52 -> 0
[ 4] 0/53 -> 0
[ 5] 0/53 -> 0
[ 6] 0/53 -> 0
[ 7] 0/53 -> 0
[ 8] 0/53 -> 0
[ 9] 0/53 -> 0
[ 10] 0/53 -> 0
[ 11] 0/53 -> 0
[ 12] 0/53 -> 0
[ 13] 0/53 -> 0
[ 14] 0/53 -> 0
[ 15] 0/53 -> 0
[ 16] 0/53 -> 0
[ 17] 0/53 -> 0
[ 18] 0/53 -> 0
[ 19] 0/53 -> 0
[ 20] 0/53 -> 0
[ 21] 0/53 -> 0
[ 22] 0/53 -> 0
[ 23] 0/53 -> 0
[ 24] 0/53 -> 0
[ 25] 0/52 -> 0
[ 26] 0/52 -> 0
[ 27] 0/52 -> 0
[ 28] 0/52 -> 0
[ 29] 0/52 -> 0
[ 30] 0/52 -> 0
[ 31] 0/52 -> 0
[ 32] 12/52 -> 0
[ 33] 40/52 -> 1
[ 34] 46/52 -> 1
[ 35] 46/52 -> 1
[ 36] 46/52 -> 1
[ 37] 46/52 -> 1
[ 38] 46/52 -> 1
[ 39] 46/52 -> 1
[ 40] 46/52 -> 1
[ 41] 46/52 -> 1
[ 42] 46/52 -> 1
[ 43] 43/52 -> 1
[ 44] 44/52 -> 1
[ 45] 35/52 -> 1
[ 46] 32/52 -> 0
[ 47] 10/52 -> 0
[ 48] 0/52 -> 0
[ 49] 0/52 -> 0
[ 50] 0/52 -> 0
[ 51] 0/52 -> 0
[ 52] 0/52 -> 0
[ 53] 0/52 -> 0
[ 54] 0/52 -> 0
[ 55] 0/52 -> 0
[ 56] 7/52 -> 0
[ 57] 42/52 -> 1
[ 58] 51/52 -> 1
[ 59] 51/52 -> 1
[ 60] 51/52 -> 1
[ 61] 51/52 -> 1
[ 62] 51/52 -> 1
[ 63] 51/52 -> 1
[ 64] 51/52 -> 1
[ 65] 51/52 -> 1
[ 66] 51/52 -> 1
[ 67] 49/52 -> 1
[ 68] 49/52 -> 1
[ 69] 42/52 -> 1
[ 70] 43/52 -> 1
[ 71] 15/52 -> 0
[ 72] 0/52 -> 0
[ 73] 0/52 -> 0
[ 74] 0/52 -> 0
[ 75] 0/52 -> 0
[ 76] 0/52 -> 0
[ 77] 0/52 -> 0
[ 78] 0/52 -> 0
[ 79] 0/52 -> 0
[ 80] 8/52 -> 0
[ 81] 47/52 -> 1
[ 82] 52/52 -> 1
[ 83] 52/52 -> 1
[ 84] 52/52 -> 1
[ 85] 52/52 -> 1
[ 86] 52/52 -> 1
[ 87] 52/52 -> 1
[ 88] 52/52 -> 1
[ 89] 52/52 -> 1
[ 90] 52/52 -> 1
[ 91] 52/52 -> 1
[ 92] 52/52 -> 1
[ 93] 39/52 -> 1
[ 94] 38/52 -> 1
[ 95] 14/52 -> 0
[ 96] 0/52 -> 0
[ 97] 0/52 -> 0
[ 98] 0/52 -> 0
[ 99] 0/52 -> 0
[100] 0/52 -> 0
[101] 0/52 -> 0
[102] 0/52 -> 0
[103] 0/52 -> 0
[104] 8/52 -> 0
[105] 44/52 -> 1
[106] 51/52 -> 1
[107] 51/52 -> 1
[108] 51/52 -> 1
[109] 51/52 -> 1
[110] 51/52 -> 1
[111] 51/52 -> 1
[112] 51/52 -> 1
[113] 51/52 -> 1
[114] 51/52 -> 1
[115] 50/52 -> 1
[116] 51/52 -> 1
[117] 43/52 -> 1
[118] 38/52 -> 1
[119] 14/52 -> 0
[120] 0/52 -> 0
[121] 0/52 -> 0
[122] 0/52 -> 0
[123] 0/52 -> 0
[124] 0/52 -> 0
[125] 0/52 -> 0
[126] 0/52 -> 0
[127] 0/52 -> 0
[128] 10/52 -> 0
[129] 46/52 -> 1
[130] 52/52 -> 1
[131] 52/52 -> 1
[132] 52/52 -> 1
[133] 52/52 -> 1
[134] 52/52 -> 1
[135] 52/52 -> 1
[136] 52/52 -> 1
[137] 52/52 -> 1
[138] 52/52 -> 1
[139] 50/52 -> 1
[140] 52/52 -> 1
[141] 42/52 -> 1
[142] 35/52 -> 1
[143] 11/52 -> 0
[144] 0/52 -> 0
[145] 0/52 -> 0
[146] 0/52 -> 0
[147] 0/52 -> 0
[148] 0/52 -> 0
[149] 0/52 -> 0
[150] 0/52 -> 0
[151] 0/52 -> 0
[152] 0/52 -> 0
[153] 34/52 -> 1
[154] 50/52 -> 1
[155] 45/52 -> 1
[156] 47/52 -> 1
[157] 43/52 -> 1
[158] 37/52 -> 1
[159] 14/52 -> 0
[160] 0/52 -> 0
[161] 0/52 -> 0
[162] 0/52 -> 0
[163] 0/52 -> 0
[164] 0/52 -> 0
[165] 0/52 -> 0
[166] 0/52 -> 0
[167] 0/52 -> 0
[168] 0/52 -> 0
[ 1] 0/52 -> 0
occInfo = cbind(uTOW,okocc)
Here are a bunch of ways to visualize the occupancy detection algorithm. The best one is all the way at the bottom, which is mostly scripted in Javascript with Charts.js.
library(ggplot2)
tow <- 40
okTOW = intervalOfWeek2==tow
ggplot() +
geom_point(data=data.frame(x=TempF,y=loadVec2), aes(x, y), color='gray') +
geom_point(data=data.frame(x=TempF[okTOW],y=loadVec2[okTOW]), aes(x, y), color='blue') +
geom_line(data=data.frame(x=TempF,y=predict(amod)), aes(x, y)) +
xlab("Temperature (°F)") + #"\xB0"
ylab("Load (kW)")
# To repeat the same data in every panel, simply construct a data frame
# that does not contain the faceting variable.
mydata=data.frame(x=TempF,y=loadVec2,ypred=predict(amod),tow=intervalOfWeek2)
mydata[mydata$tow<10,]
ggplot(mydata[mydata$tow<10,], aes(x, y)) +
geom_point(data = transform(mydata, tow = NULL), colour = "grey85") +
geom_point() +
facet_wrap(vars(tow))
library(plotly)
Attaching package: 㤼㸱plotly㤼㸲
The following object is masked from 㤼㸱package:ggplot2㤼㸲:
last_plot
The following object is masked from 㤼㸱package:stats㤼㸲:
filter
The following object is masked from 㤼㸱package:graphics㤼㸲:
layout
mydata %>%
plot_ly(
x = ~x,
y = ~y,
frame = ~tow,
type = 'scatter',
mode = 'markers',
showlegend = F
)
# geom_point(aes(x, y), color='gray') +
# geom_line(aes(x, ypred)) +
p <- ggplot(mydata) +
geom_point(aes(x, y, frame=tow), color='blue') +
xlab("Temperature (°F)") + #"\xB0"
ylab("Load (kW)")
Ignoring unknown aesthetics: frame
fig <- ggplotly(p) %>%
animation_opts(
frame=500, easing='linear', redraw = F
)
fig
This suffers from slow speed (heavy on the CPU usage).
mydata %>%
plot_ly(x=~x,mode='none') %>%
layout(xaxis=list(title=list(text='Temperature (°F)')),
yaxis=list(title=list(text='Load (kW)'))) %>%
add_trace(y=~y,name='loads', type='scatter', mode='markers',
showlegend=F, marker=list(color='grey'), text = 'none') %>%
add_lines(y=~ypred,name='fit',
showlegend=F, line=list(color = 'black')) %>%
add_trace(
y = ~y,
frame = ~tow,
type = 'scatter',
mode = 'markers',
showlegend = F,
marker = list(color = 'blue')
) %>%
animation_opts(
frame=500, easing='linear', redraw=F
)
Next, I’ve created a javascript-based charting widget because this is more interactive than an animation, and works better because it doesn’t animate the point cloud in the background. The widget does not yet interact with R, so I’m just going to show how to export some data to a file that I have manually pasted into the widget. Here are some future features:
Here is how I export the data to the widget.
# Convert to JSON for Chart.js
dataByTOW=uTOW %>% purrr::map(~ mydata[mydata$tow==.x,])
dataByTOWjson = jsonlite::toJSON(dataByTOW, pretty=F)
#fitLine = unique(mydata[with(mydata,order(x)),c('x','ypred')])
fitLine = data.frame(x=TempF,y=predict(amod)) %>% unique() %>% arrange(x)
fitLinejson = jsonlite::toJSON(fitLine, pretty=F)
sink('datafile.js')
cat(glue('const dataByTOW = {dataByTOWjson};
const fitLine = {fitLinejson};'))
sink()
#mydata %>% mutate(x=x,y=ypred)
Here is the widget. Try pointing the mouse over the time-of-week grid:
That is all. Thank you for reading. Find more from my fork of RMV2.0